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The evolution of AGN across cosmic time: what is downsizing? 
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ABSTRACT 

We use a coupled model of the formation and evolution of galaxies and black holes (BH) 
to study the evolution of active galactic nuclei (AGN) in a cold dark matter universe. The 
model predicts the BH mass, spin and mass accretion history. BH mass grows via accretion 
triggered by discs becoming dynamically unstable or galaxy mergers (called the starburst 
mode) and accretion from quasi-hydrostatic hot gas haloes (called the hot-halo mode). By 
taking into account AGN obscuration, we obtain a very good fit to the observed luminosity 
functions (LF) of AGN (optical, soft and hard X-ray, and bolometric) for a wide range of 
redshifts (0 < z < 6). The model predicts a hierarchical build up of BH mass, with the typical 
mass of actively growing BHs increasing with decreasing redshift. Remarkably, despite this, 
we find downsizing in the AGN population, in terms of the differential growth with redshift 
of the space density of faint and bright AGN. This arises naturally from the interplay between 
the starburst and hot-halo accretion modes. The faint end of the LF is dominated by massive 
BHs experiencing quiescent accretion via a thick disc, primarily during the hot-halo mode. 
The bright end of the LF, on the other hand, is dominated by AGN which host BHs accreting 
close to or in excess of the Eddington limit during the starburst mode. The model predicts 
that the comoving space density of AGN peaks at z ~ 3, similar to the star formation history. 
However, when taking into account obscuration, the space density of faint AGN peaks at lower 
redshift (z < 2) than that of bright AGN (z ~ 2 — 3). This implies that the cosmic evolution 
of AGN is shaped in part by obscuration. 
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1 INTRODUCTION 



Reproducing the evolution of active galactic nuclei (AGN) is 
a crucial test of galaxy formation models. Not only are AGN 
an important component of the Universe, they are also thought 
to play a key role in the evolution of their host galaxies and 
their environment. For example, feedback processes which ac- 
company the growth of black holes (BHs) during AGN ac- 
tivity are thought to influence star formation (SF) through 
the suppression of cooling flows in galaxy groups and clus- 
ters (Dalla Vecchia et al. 2004; Springel, Di Matteo, & Hernquist 
20051- ICroton et al.ll2006l; iBower et alj|2006l: iHopkins et alj|200a 
Thacker, Scannapieco, & Couchman 2006; Some rville et al.l 



standing these processes requires that the underlying galaxy forma- 
tion model is able to track and reproduce the evolution of AGN. 



Early surveys demonstrated that 



quasi-stellar objects 
3m z ~ up to 



- 2.5 dSchmidt & Green 


1 19831: iBovle, Shanks & Peterson] 


Hewett, Foltz. & Chaffee 


19931; Bovleetal. 2OO0l). Be- 



yond z 



2, the space density of QSOs starts to decline 



( Warren. Hewett, & Osmerl 1 1 994 


; |Schmidt, Schneider, & Gunnl 


1 19951; Fan et al. 200 ll; IWolf et al. 


2003). Recently there has been 



Lagos. Cora, & Padilld 120081 see also iMarulli et al.l l2008h . Sim 



ilarly, energy feedback from AGN is a key ingredient in 
determining the X-ray properties of the intracluster medium 
(Bower, Mc Carthy. & Benson 2008; McCart hy et alj20ld) . Under- 



much progress in pinning down the evolution of faint AGN in 
the optical. Using the 2-degree-field (2dF) QSO Redshift survey, 
the luminosity function (LF) was prob ed to around a magnitud e 
fainter than the break to z ~ 2 (2QZ lCroom et al.ll200lU2004h . 
This limit was extended by t he 2dF-SDSS lumino u s red galaxy 
and Q SO (2SLAQ) survey dRichards et all 120051 ; ICroom et all 
2009a). The estimate of the LF fr om the final 2SLAQ catalogue 
reached M b ~ -19.8 a t z = .4 dCroom et al]|2009bh . With the 
2SLAQ LF. ICroom et alj j2009bl) demonstrated that in the optical, 
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faint quasars undergo mild evolution, with their number density 
pea king at lower redshifts than is the case for bright quasars (see 

also lBongiorno et al.ll2007T) . 

F aint AGN can be sele c ted robustly in X-ray s (Hasing er et al 



2001 



iGiacconi et al.| [2002; Ale xander et al J [20031; iBarcons et al 



20071 see also the review by Brandt & Hasinger 2005), which 
means that a wider range of AGN luminosity can be probed in X- 
rays than is possible in the optical. This permits the study of the 
evolution of a wider variety of AGN in addition to QSOs (e.g., 
Seyfert galaxies) and thus provides a more representative picture 
of the various AGN populations. The evolution of the AGN LF in 
X-rays has been investigated by many authors by employing data 
from th e Chandra, ASCA, ROSAT, HEAO-1 and XMM-Newton 
surveys |Mivaii. Hasinger. & Schmidtl2000l:lLa Franca et al .120021; 



Cowie et al.ll2003l;lFiore et alJl2003|;IUeda et alJl2003l ; lBarger et all 
2005; Hasinger. Mivaii. & Schmidtll2005l) . These studies show that 
in both soft (0.5 - 2keV) and hard (2 - lOkeV) X-rays, faint AGN 
are found to evolve modestly with redshift. In contrast, bright AGN 
show strong evolution, similar to that seen for quasars in the optical. 
In addition, observations in soft X-rays suggest that the comoving 
space density of brigh t AGN peaks at higher redshifts (z ~ 2) than 
faint AGN (z < l. lHasinger. Mivaii. & Schmidtll2005h . 

The differential evolution of bright an d faint AGN with 
redshift has been described as d ownsizing dBarger et al.l 120051 : 
lHasinger. Mivaii. & Schmidtl 120051) . This implies that AGN ac- 
tivity in the low-z universe is dominated by either high-mass 
BHs accreting at low rates or low-mass BHs growing rapidly. 
Hopki mTet al.l l|2005b) proposed that the faint end of the LF is 
co mposed of high mass B H s experiencing qui escent accretion (see 
also Hopkins et al.||2005 JcllBabic et alj|2007l) . The bright end, on 
the other hand, in this pictur e corresponds to BHs accreting near 
their Eddington limit. In the iHopkins et alj model, quasar activ- 
ity is short-lived and is assumed to be driven by galaxy mergers 
dPi Matteo. Springel. & He rnquist 200^). 

The mass of accreting BHs can be estimated using the spectra 
of the AGN. Quasars are ideal for this since they are detected up 
to very high redshifts and hence can trace the evolution of actively 
growing BH mass back into the early Universe (Fan et al. 2001; 
Fan et al. 2003; Fontanot et al. 2007; Willott et al. 2010) . The 
BH mass in quasars is calculated using empirical relations derived 
from optical or UV spectroscopy. In particular, mass-scaling re- 
lations between the widths of different broad emission lines and 
continuum luminosities that have been ca l ibrated against reverber- 
ation mapping results dVestergaardl |2002| ; IVestergaard & Petersonl 
120061) have allowed BH mass estimates in several large, unobscured 
(type-1) quasar samples. When translating quasar luminosities into 
BH masses using t he width of broad lines, a similar downsizing i s 
seen in BH mass dVestergaard & Osmed[2009l ; iKellv et alJboiOh . 
suggesting that the most massive BHs (A/bh > 10 9 Mq) were al- 
ready in place at z > 2, whereas the growth of the less massive 
ones is delayed to lower redshifts. 

An important feature of AGN which should be taken into ac- 
count is that they exhibit evidence of obscuration at both optical and 
soft X-ray wavelengths. The obscuration may be linked to the exis- 
tence of a geometrical torus around the accretion disc whose pres- 
ence is invoked in the AGN unification scheme jAntonucci|[l993l ; 
lUrrv & Padovanil 1 1995b or to intervenin g dust clouds related to 
physi cal processes within th e host galaxy jMartmez-Sansigre et al.l 
2005), s uch as SF activity feallantyne. Everett. & Murray! 120061 
but also iGoulding & Alexanderj|2009l) . As a consequence, a large 
fraction of AGN could be obscured and thus missing from opti- 
cal and soft X-ray surveys. Hence, when applying the scaling re- 



lations to estimate the mass of accreting BHs, the absence of ob- 
scured (type-2) quasars from the AGN samples may introduce sig- 
nificant biases into the inferred evolution of BH mass. Only hard 
X-rays can directly probe the central engine by penetrating the ob- 
scuring medi um, therefore providing complete and unbiased sam- 
ples o f AGN dUeda et al.ll2003l ; lLa Franca et"ai]|2005l ; iBarger et al.l 
2005). However, even in hard X-rays, a population of Compton- 
thick sources, namely AGN with hydrogen colu mn densities ex- 
ceeding Nh ~ 10 cm" , would still be missing (Comastri2004; 
lAlexander et alj2005ll2008l ; lGoulding et alfcOld) . 

In this paper we pr esent a study o f AGN evolution using semi- 
analytic modelling (see Baugh] |2006l for a review). Our aim is to 
provide a robust framework for understanding the downsizing of 
AGN within a self-consistent galaxy formation model. The paper 
is organized as follows. In Section 2 we present the galaxy forma- 
tion model upon which we build our AGN model. In Section 3 we 
study the evolution of BH mass and explore the scaling relations 
predicted between the BH mass and the mass of the host galaxy 
mass and dark-matter (DM) halo. In Section 4 we present the es- 
sential ingredients of the AGN model and study the evolution of 
the physical properties which, together with the mass, provide a 
complete description of accreting BHs. In Section 5 we present 
our predictions for the evolution of the optical, soft/hard X-ray and 
bolometric LFs. Finally, in Section 6 we explore the downsizing of 
AGN in our model. 



2 THE GALAXY FORMATION MODEL 

We use the GALFORM se mi-analytical galax y formation code 



I galaxy 

dCole etalj|l994ll2000l. also lBaugh et ai]|2005l ; lBower et alj|2006l; 
iFont et all 20081) and the extension to follow the evolution of BH 
mass and spin introduced by Fanidakis et al. (2010). GALFORM 
simulates the formation and evolution of galaxies and BHs in a 
hierarchical cosmology by modelling a wide range of physical pro- 
cesses, including gas cooling, AGN heating, SF and supernovae 
(SN) feedback, chemical evolu tion, and galaxy me rgers. 

Our starting point is the iBower et al.l d2006h galaxy forma- 
tion model. This model invokes AGN feedback to suppress the 
cooling of gas in DM haloes with quasi-static hot atmospheres 
and has been shown to reproduce many observables, such as 
galaxy colours, stellar masses and LFs remarkably well. The 
model adopts a BH growth recipe based on that introduced by 
Malb on et al] d2007l) . and extended by Bower et al. (2006) and 
Fanidaki s et al. (2010). During starburst s triggered by a disc in- 
stability (Efstath iou. Lake & NegroponteHl982r) or galaxy merger, 
the BH accretes a fixed fraction, /bh, of the cold gas that is turned 
into stars in the burst, after taking into account SN feedback and 
recycling. The value of /bh is chosen to fit the amplitude of the 



channel, quiescent inflows of gas from the hot ha lo during AGN 
feedback also contribute to the mass of the BH (see lWhite & Frenkl 
ll99ll ; ICole e"tai]|2000l; ICroton et alj|200d for the cooling proper- 
ties of gas in haloes). To distinguish between these two accretion 
channels, we refer to the accretion triggered by disc instabilities or 
galaxy mergers as the starburst mode and the accretion from quasi- 
hydrostatic hot haloes as the hot-halo mode. These correspond to 
th e quasar and r adio m odes, respectively, in the terminology used 
bv lCrotonetai1d2006l) . Finally, mergers between BHs, which oc- 
cur when the galaxies which host the BHs merge, redistribute BH 
mass and contribute to the build up of the most massive BHs in the 
universe (see Fanidakis et al. 2010). 
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Table 1. Summary of the revised parameter values in the variants of the 
Bower et al. model considered here. 



Study 


r a 
JEdd 


2 b 
e kin 


/bh c 


Fan 10b (this study) 


0.039 


0.016 


0.005 


FanlOa (Fanidakis et al. 2010) 


0.01 


0.1 


0.017 



Notes. "Fraction of the Eddington luminosity available for heating the gas 
in the host halo. 6 Average kinetic efficiency during the hot-halo mode. 
c Ratio of the mass of cold gas accreted onto the BH to the mass of cold 
gas used to form stars. 



GALFORM calculates a plethora of properties for each galaxy 
including disc and bulge sizes, luminosities, colours and metallici- 
ties to list but a few. In this analysis, we are primarily interested in 
the properties that describe the BHs. Specifically, we output the BH 
mass, Mbh, the amount of gas accreted in a starburst or during the 
hot-halo mode, Af acc , the time that has passed since the last burst 
of SF experienced by the host galaxy and the BH spin (discussed 
in more detail below). The time since the last burst is necessary to 
determine the beginning of the active phase of BH growth in the 
starburst mode. 

Fanidakis et al. (2010) extended the model to track 
the evolution of BH spin, a. By predicting the spin we 
can compute interesting properties such as the efficiency, e, 
of converting matter into radiation during the accretion of 
gas onto a BH and the me chanical energy of jets in the 
iBlandford & Znaiekl d 19771) and lBlandford & Pavnd d 19821) mech- 
anisms. Toj^ri£yilateThe_em use the chaotic accretion 
model llKing. Pringle. & Hofmannl 120081 : iBerti & Volonterill2008l : 
iFanidakis et alT feoiO). In this model, during the active phase gas 
is accreted onto the BH via a series of randomly orie nted accre- 
tion d iscs, whose mass is limited by their self-gravity dKing et al.l 
l2005h . We choose the cha otic model (over the alternative prol onged 
model in the literature bv lVolonteri. Sikora. & La sota 2007) since 
in this case the predicted BH spin distributions, through their influ- 
ence on the strength of relativistic jets, reproduce better the pop- 
ulatio n of radio-loud AGN in the local Universe (Fanida kis et al.l 
l20ld) . A correction to the amount of gas that is accreted is ap- 
plied in order to account for the fraction of gas that turns into radi- 
ation during the accretion process. Finally, we note that we do not 
take into accoun t BH ejection via gravitational-wave recoils during 
galaxy mergers dMerritt et alj2004h . This omission is not expected 
to hav e a significant impact on the evolution of BH mass in our 
model ( Libeskind et al. 20061). 

With this extension to the iBower et alj model, in 
IFanidakis et alj we were able to reproduce the diversity of 
nuclear activity seen in the local universe. Furthermore, we 
demonstrated that the bulk of the phenomenology of AGN can be 
naturally explained in a ACDM universe by the coeval evolution of 
galaxies and BHs, coupled by AGN feedback. In this paper we aim 
to extend the predictive power of the model to high redshifts and 
different wavelengths, to provide a complete and self-consistent 
framework for the formation and evolution of AGN in a ACDM 
cosmology. 

For completeness, we now list the model parameters which 
have an influence over the growth of BH mass: 

(i) /bh, the fraction of the mass of stars produced in a star- 
burst, after taking into account SN feedback and the recycling of 
gas, that is accreted onto the BH during the burst. 

(ii) fvdd, the fraction of the Eddington luminosity of an ac- 




Figure 1. Top: The cosmic history of the total SFR density (black lines) in 
the FanlOb (solid line) and FanlOa (dashed line) models. Also shown is the 
contribution from each SF activity mode, namely burst and quiescent SF 
(red and blue respectively), to the total SFR density. Bottom: The contri- 
bution of disc instabilities (green lines) and galaxy mergers (orange lines) 
to the cosmic history of SFR density in the burst SF activity mode for the 
FanlOb (solid lines) and FanlOa (dashed lines) models. 



creting BH that is available to heat the hot halo during an episode 
of AGN feedback. 

(iii) 6ki n , the average kinetic efficiency of the jet during the 
hot-halo mode. 

The fiducial model in this paper is denoted as FanlOb. We 
adopt the values for the above parameters that were used in the 
Bower et al. (2006) mo del, as listed in Table 1 against the FanlOb 
model. Fanidaki s et alj made some small changes to these param- 
eter values, which are given in the FanlOa entry in Table 1. These 
changes resulted in a small change in the distribution of accretion 
rates in the hot-halo mode. By reverting to the original parameter 
values used in the Bower et al. model, we note that there is a small 
tail of objects accreting in the hot-halo mode, for which the ac- 
cretion rate is higher than that typically associated with advection 
dominated accretion via a thick disc (see Section 4.1 for further 
discussion). 

A further difference between our fiducial model, FanlOb, and 
the FanlOa model is the use of an improved SF law. Following 
Lagos et al. (2010), who implemented a range of empirical and 
theoretical SF laws into GALFORM, we use the SF law of Blitz & 
Rosolowsky (2006, hereafter BR06). Lagos et al. found that this 
SF law in particular improved the agreement between the model 
predictions and the observations for the mass function of atomic 
hydrogen in galaxies and the hydrogen mass to luminosity ratio. 
The BR06 model has the attraction th at it is more physical than the 
previous parametric SF law used in lBower et all and agrees with 
observations of the surface density of gas and SF in galaxies. The 
BR06 law distinguishes between molecular and atomic hydrogen, 
with only the molecular hydrogen taking part in SF. The fraction 
of hydrogen in molecular form depends on the pressure within the 
galactic disc, which in turn is derived from the mass of gas and 
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z 

Figure 2. The fraction of baryons locked in BHs as a function of redshift 
(black lines) for the FanlOb (solid lines) and FanlOa (dashed lines) models. 
The plot also shows the contribution of each accretion channel, namely disc 
instabilities (orange), galaxy mergers (green) and quasi-hydrostatic halo ac- 
cretion (blue), to the total fraction of baryons in BHs. 

stars and the radius of the disc; these quantities are predicted by 
GALFORM. The BR06 SF law contains no free parameters once it 
has been calibrated against observations. We adopt a value for the 
normalization of the surface density of SF that is two thirds of the 
value used by Lagos et al., but which is still within the observa- 
tional uncertainty. We do this to improve the match to the observed 
quasar luminosity function. 

As shown by Lagos et al., the adoption of the BR06 SF law 
changes the SF history predicted by the model. This is due to a 
much weaker dependence of the effective SF timescale on redshift 
with the new SF law, compared with that displayed by the Bower 
et al. model. This leads to the build-up of larger gas reservoirs in 
discs at high redshift, resulting in more SF in bursts. This is shown 
in the top panel of Fig. Q] which compares the predictions of the 
FanlOb (solid lines) and FanlOa (dashed lines) models for the to- 
tal SFR density and the SFR density in the quiescent and burst SF 
modes. Fig.Q]shows that the individual SF modes change substan- 
tially on using the BR06 SF law. The quiescent SF mode in the 
FanlOb model is suppressed by almost an order of magnitude above 
z ~ 3 compared to the FanlOa model. This has an impact on the 
build up of BH mass, by changing the amount of mass brought in 
through the starburst mode. The enhancement of the burst mode is 
further demonstrated in the bottom panel of Fig.Q] where we show 
the SFR density history in bursts distinguishing between those trig- 
gered by disc instabilities and galaxy mergers for the FanlOa and 
FanlOb models. Both channels show a significant increase in SFR 
density in the FanlOb model. Since the BHs in our model grow in 
bursts of SF following disc instabilities and galaxy mergers we ex- 
pect this enhancement to have a significant impact on the evolution 
of BH mass. 

To evaluate the change in the BH mass when we use the BR06 
SF law, we plot in Fig.[2]the fraction of baryons locked up in BHs 
as a function of redshift. In this plot we have also distinguished 



Table 2. Summary of the local BH mass densities found in this and previous 
studies (assuming h = 0.7). Densities are shown for BHs in disc (S,S0) and 
elliptical (E) galaxies, and for the global BH population in our model (tot). 



Study pbh^SO) 1 pbh(E) 2 p BH (tot) 3 



This study 


0.72 


4.70 


5 


Graham et al. (2007) 


95+ ' 49 
u - yo -0.49 


3 46+ 1 16 


4.41 


Shankar et al. (2004) 


i 1+0.5 

^-o.s 


q 1+0.9 

■^-o.s 


4.2 


Marconi et al. (2004) 


1.3 


3.3 


1.6 


Fukugita & Peebles (2004) 


1 7+ 1 ' 7 
1 -'-0.8 


3.41?;* 


5.1 



Notes. 1 ' 2 > 3 In units of 10 s M Q Mpc~ 3 . 



between the different accretion channels: disc instabilities, galaxy 
mergers and accretion from hot gas haloes in quasi-hydrostatic 
equilibrium. We note that in this plot it becomes clear that the 
growth of BHs is dominated by accretion of cold gas during disc 
instabilities. Therefore, secular processes in galaxies are responsi- 
ble for building most of the BH mass, while galaxy mergers be- 
come an important channel only when they occur between galaxies 
of similar mass (i.e. major mergers). 

The contribution of the starburst mode to the build up of BH 
mass increases significantly in the FanlOb model compared with 
the FanlOa model. Both the disc instability and galaxy merger 
channels are enhanced by a significant factor, resulting in differ- 
ent evolution of the total BH mass. The FanlOb model also predicts 
different hot-halo mode accretion. This can be attributed to the dif- 
ferent values of the ekin and /Edd parameters used in FanlOa (note 
that halo accretion is completely independent of SF in the galac- 
tic discs). It will become evident in later sections that all of the 
aforementioned changes induced to the BH mass are essential for 
providing a good fit to the evolution of the AGN LF. 

Finally we note that the total SFR density history in Fig. Q] 
remains unchanged when using the BR06 SF law, as reported by 
Lagos et al. This implies that properties such as stellar masses and 
galaxy luminosities are fairly insensitive to the choice of SF law. 
Indeed, lLagos et"al] show that the predictions for the &j and K- 
band LF s on using the BR06 SF law are very similar to those of the 
original iBower et al.l model and reproduce the observations very 
closely. Since we have changed the normalization of the SF surface 
density, we recheck that our model still reproduces the observed LF 
of local galaxies, as shown in Fig. [3] Hence, we are confident that 
we are a building an AGN model within a realistic galaxy formation 
model. 



3 THE EVOLUTION OF BH MASS 

BH demographics have been the topic of many studies in the past 
decade mainly because of the tight correlations between the prop- 
erties of BHs and their host stellar spheroids. These correlations 
take various forms, relating, for example, the mass of the BH 
to the mass of the ga l actic bulge (the Mbh — MbuIrb relation: 
i Magprrian et alj|l998t iMcLure & Dunlodl2002l ; iMarconi & Hunt! 
booilHaring & Rixll2004l), or to the stellar velocity dispersion (the 
Mbh - o- relation: iFerrarese & Merrittl200d : lGebhardt et alj2OO0l: 
Trem aine et "al] |2002h . These remarkable and unexpected correla- 
tions suggest a natural link between the BH evolution and the for- 
mation history of galaxies. The manifestation of this link could be 
associated with AGN activity triggered during the build up of BHs. 
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Figure 3. The LF of galaxies in the local Universe. The left panel compares the predictions of the FanlOb (solid line) and FanlOa (dashed line) models for the 
6j-band LF with the observational determination from the 2dF galaxy redshift survey by Norberg et al. (2002). Similarly, the right panel shows the predictions 
of the models (line style as before) for the K-band LF and compared to the observational determinations by Cole et al. (2001) and Kochanek et al. (2001). The 
theoretical predictions from both models include dust extinction. 



The observed Mbh — MB U i gc relation can be used to esti- 
mate the mass of a BH. Several authors have utilised this technique 
to estimate the BH m a ss function (BHMF) in the local Univers e 
dYu&Tremaiiiell2002l ; iMarconi et al]|2004 IShankar et ai]l2004h . 
using scaling relations such as A4bh — o", and A/bh — i/Buige- 
Based on these MFs, and others inferred by independent studies, 
the total BH mass density in the local Universe has been estimated 
to be in the range p B H = (4.2 - 5.1) x 10 5 M Mpc" 3 (see 
Table [2] for a list of local BH mass density estimates). Note that, 
when correcting for unknown dependencies on h, the different esti- 
mates from the literature sugg est a BH mass densit y consistent with 
4.4 - 5.9 x 10 5 M Mpc" 3 jGraham et ai1l2007l) . 

The global BHMF predicted by the model for all BHs at z = 
is shown in Fig.|4^. The predicted MF is almost constant in the mass 
range 10 6 — 10 s M@. For higher masses the MF decreases steeply 
with increasing mass. Our predictions are compared to the loca l 
MF estimated by IMarconi etal] d2004) and IShankar et al] J2004l) 
using the BH mass scaling relations in local galaxies. The pre- 
dicted and observed MFs are in good agreement for low mass BHs 
and disagrees moderately for the 10 9 - 10 10 M BHs. The over- 
all BH mass density at z = predicted by our model amounts to 
p BH = 5.42 x 10 5 M Mpc -3 and is consistent with that estimated 
bv lMarconi et alJfflRH = 4.6li ! x 10 5 M B Mpc~ 3 ) and other au- 
thors JShankar etal]|2004l : lGraham et alJl2007l : lYu & Lull2008l, see 
Table O. Note that, when calculating the BH mass density we do 
not take into account mas s loses due to gravitational wave emission 
dMenou & Haimanll2004h . 

In Fig.|4^ we also show the contribution to the MF from BHs 
hosted by elliptical galaxies. We classify as ellipticals all the galax- 
ies whose bulge lumin osity contributes more than 60% to their total 
6,j-band luminosity JCole et al. 2000; Parry et al. 2009). The rest of 
the galaxies are classified as spiral or SO galaxies. As expected, the 
BHs in the elliptical galaxies dominate the high mass end of the MF. 
Hence, the most massive BHs inhabit the centres of the most mas- 
sive systems in our simulations. By contrast, spiral and SO galaxies 
contribute only to the low-mass end of the MF. 



The growth of BHs in spiral and elliptical galaxies is driven in 
principle by different channels. In spiral galaxies, BHs grow mainly 
via accretion during the starburst mode whereas the massive BHs in 
ellipticals grow via accretion during the hot-halo mode or BH-BH 
mergers. Evidently, the evolution of the BHMF should be shaped 
by these growth channels. This is illustrated in Fig.[4jj, where we 
show the evolution of the global BHMF in the redshift range z = 
— 6.2. The build up of the BHs populating the low-mass end takes 
place at high redshifts and is complete at 2 ~ 2. Since z ~ 2 the 
amplitude of the BHMF remains almost unchanged indicating that 
10 6 — 10 s Mq BHs in spiral galaxies today were already in place 
by z ~ 2, which follows from the spheroid being older than the 
disc. The growth of these BHs is dominated by accretion during 
the starburst mode. 

The fact that the space density of this population does not 
evolve since z ~ 2 suggests that the starburst mode becomes less 
significant in the low-redshift universe. Nonetheless, the build up 
of the BH mass continues below z ~ 2, where we witness the 
build up of the most massive BHs in the universe. These BHs grow 
mainl y during the hot-hal o mode or via mergers with other BHs 
(see iFanidakis et al]|2010l) . The different physics governing these 
channels results in a different slope for the high-mass end of the 
BHMF giving rise to a strong break at ~ 5 x 10 s Mq. Yet, the 
build up of these BHs is not very efficient. This becomes clearer 
when we consider the evolution of the BH mass density. For in- 
stance, in the redshift range z — 6.2 — 1.9 the BH mass density 
increases by a factor of ~ 25, as a result of the intense accretion 
during the starburst mode in the high redshift universe. By contrast, 
from z — 1.9 until z = the density only doubles, mainly because 
the 10 9 Mq BHs become more numerous. 

To gain more insight into the rapidly evolving population of 
accreting BHs we consider the evolution of actively growing BHs, 
namely those BHs that accrete at relatively high rates (greater that 
0.1% of their Eddington accretion rate, see Section|4j- In Fig.|4j; 
we show the MF of actively growing BHs in our model and its 
evolution with redshift. The most striking characteristic of the MF 
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Figure 4. a) The global BHM F at z = pre dicted by our model (black solid line). The observationally determined MFs are taken from Marconi et al. (2004, 
red filled triangles) and Shankar et al. 12004, blue filled circles). The error bars correspond to ±lcr uncertainties, b) The evolution of the global BHMF with 
redshift. c) The MF of actively growing BHs (BH accreting in the thin-disc regime) calculated at different redshifts as indicated in the legend. 



is the dramatic change in shape and normalization with redshift. 
At z — the MF has a maximum at ~ 10 7 Mq below which 
the accreting BHs have an almost cons tant space density (see also 
results from SDSS lHeckman et alj|2004h . A significant bend in the 
slope located at ~ 10 8 Mq is seen, a feature that remains evident 
up to z ~ 2 — 3. At higher redshifts, the space density of BHs 
drops substantially, establishing that accreting BHs with masses > 
10 9 Mq are very rare objects. These BHs are accreting at relatively 
low rates and therefore power faint- AGN activity (see discussion in 
Sectionl44t. 

As illustrated by Fig. [4}:, the low-mass end of the MF evolves 
differently with redshift compared to the high-mass end. Its am- 
plitude increases significantly with increasing redshift, suggesting 
that the space density of low-mass accreting BHs was higher at ear- 
lier epochs. This evolution is consistent up to z ~ 2.5 above which 
the amplitude of the low-mass end starts to decline modestly. The 
high-mass end shows a similar increase towards higher amplitudes 
with increasing redshift, however, the decline appears earlier, al- 
most at z ~ 1, This is because the high-mass end forms much later 
than the low-mass end in our model. 

The change of the BHMF with redshift illustrated in Fig. |4j; 
indicates an evolutionary scenario for the actively growing BH pop- 
ulation in our model which mirrors the hierarchical growth of their 
host galaxies. At high redshifts the accretion activity is dominated 
by 10 6 Mq BHs. Accretion onto these BHs results in the fast build 
up of the 10 7 — 10 8 Mq BHsatz ~ 4—6, as implied by the increase 
in the space density of these BHs in Fig. [4};. From z ~ 4 the accre- 
tion activity, gradually shifts to the 10 7 — 10 s Mq BH population. 
This is mainly because more higher mass galactic systems are now 
in place and thus disc instabilities and galaxy mergers contribute to 
the growth of higher mass BHs compared to earlier epochs. Accre- 
tion onto these BHs leads to the fast build up of the > 10 s M Q BH 
population, which will later be promoted via accretion during the 
hot-halo mode and BH-BH mergers to the most massive BHs in our 
model. 

Eventually, the different physical processes that drive the evo- 
lution of BH mass result in a tight correlation between the BH and 
host galaxy mass, as shown in Fig. [5^. In this figure, we plot the 
median of the log Mbh — log Msuige distribution (solid lines) and 



the associated 10—90 percentiles (dotted lines) at z = 0, 1.08, 2.42 
and 6.2. Surprisingly, a well defined Mbh — MB U i ge relation is es- 
tablished already at z = 6.2, merely as a consequence correlation 
of the accretion process with the SF that shapes the mass of the 
host bulges. Even when the growth of BHs is dominated by other 
processes, which are only indirectly linked to the SF, the tight cor- 
relation remains. For example, the most massive BHs correlate with 
the mass of their host bulge because feedback energy when gas is 
accreted onto them regulates the cooling flows that supply gas for 
SF and therefore the change of mass of the bulge. 

In a hierarchical universe the most massive galactic bulges are 
usually found in the most massive DM haloes. Since the galaxies 
with the most massive bulges (namely the passive ellipticals) host 
the most massive BHs and those with less massive bulges (namely 
the spiral and SO galaxies) host the low-mass BHs, it is natural to 
assume that there must be a similar hierarchy in the halo environ- 
ments where these BHs can be found. Indeed, when plotting the 
relation between the BH mass and the mass of the host halo we 
find a tight correlation. This is illustrated in Fig. [5]) where we show 
the median of the log Mnaio — log Mbh distribution (solid lines) 
and its 10 — 90 percentiles (dotted lines) for different redshifts. Re- 
markably, the physical processes that shape the Mauigc — Mbh 
relation in our model give rise to a well defined link between these 
two quantities at all redshifts. The A/naio — Mbh relation can be 
parametrized by a power law with a very steep slope that becomes 
shallower above M Ha io ^ 10 12 - 10 13 M s . 

The different slopes indicate the different efficiency with 
which BHs grow in haloes of different mass. Low mass haloes 
are very efficient in growing BHs since in these environments 
the gas cooling is not suppressed. As a consequence, BHs dou- 
ble their mass remarkably quickly resulting in a steep slope for the 
log Afnaio — log Mbh relation. The fast BH mass build up slows 
down when hot gas in the host haloes enters the quasi-hydrostatic 
regime. In this case, AGN feedback suppresses the cooling flows 
that provide fresh cold gas to the galaxies and thus accretion dur- 
ing the starburst mode is reduced. In the quasi-hydrostatic regime 
gas accretion during the hot-halo mode and mergers dominate the 
BH mass build up. The hot-halo mode is, however, characterized 
by low accretion rates (see Section l4~3b . BH mergers rather than 
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Figure 5. The median of the A-fBuigc — A/bh ( a ) ar >d Afjjalo — A-^bh (b) distributions predicted by our model (solid lines). Predictions are shown for 
z = 0, 1.08, 2.42 and 6.20 as indicated by the legend. The dotted lines indicate the 10 — 90 percentile ranges of the distributions. 



adding new baryons to the BHs, only redistribute the BH mass. 
Therefore, the BH mass build up slows down establishing haloes of 
A^Haio ~ 10 13 — 10 15 M Q as environments where BH growth is 
not very efficient. 

Given the efficiencies that characterize the two different 
regimes in Fig. |5p we expect to find the brightest AGN in the 
10 11 — 10 13 Mq haloes. This perhaps is in contrast with the com- 
mon expectation that the brightest quasars should be found in the 
most massive haloes. 



to the dynamical timescale of the host bulge, £Bui gc , through the 
relation 



— /q^Bu 



lge- 



(2) 



/ q is a free parameter and its value is fixed to 10. This is determined 
by fitting the predictions of our model for the quasar LF to the 
observations (see Section [5^2l (. 

A second important quantity associated with every BH is the 
Eddington luminosity, 



4 BH SPINS, ACCRETION EFFICIENCIES AND DISC 
LUMINOSITIES 

AGN are unambiguously the most powerful astrophysical sources 
in the Universe. The ability of an accretion disc to produce elec- 
tromagnetic radiation is attributed to gravity yet its radiative effi- 
ciency is controlled primarily by the properties of the gas. When 
the gas settles itself into a thin, cool, optically-thick accretion disc 
( Shak ura & Sunvaevl 19731) . the efficiency of converting matter into 
radiation can reach 32% (Novikov & Thome 1973). When the flow 
is characterized by very high, Eddington accretion rates, such ex- 
tremely luminous discs can power 10 47 — 10 48 erg s _1 quasars. In 
the next sections we describe how we model the physics of accre- 
tion flows and present predictions for the most fundamental prop- 
erties characterising the accreting BH systems in our model. 



4.1 Calculation of the disc luminosity 

The first important property that we can calculate in our model is 
the physical accretion rate, M, onto a BH. This is defined as 

A^/acc 



M 



t. 



(1) 



where Af aC c is the total accreted mass and t acc is the accretion 
timescale. The accretion timescale is assumed to be directly linked 



4-KGM B nm p c 46 

-t^Edd = = 1.4 X 1U 



Mbh 
10 8 M 



erg s 



(3) 

where k ~ 0.3 cm 2 g _1 is the electron scattering opacity and m p 
is the proton mass. The Eddington luminosity has an associated 
accretion rate which is expressed as 



A/sdd = L Edd /ec 2 



(4) 



where e is the accretion efficiency and c is the speed of light. This is 
the accretion rate at which the black hole radiates at the Eddington 
luminosity. The accretion efficiency is assumed to be determined 
by the spin of the BH (Novikov & T horne 1973) and its value is 
calculated as in iFamdalris et ailboid) . It is convenient in our anal- 
ysis to express the physical accretion rate in units of the Eddington 
accretion rate, m = M /M^dd, in order to introduce a dependence 
on the BH mass. 

The accretion rate has a dramatic impact on the ge- 
ometry and radiative propertie s of the accretion disc (see 
iDone. Gierlinski. & Ku bota 2007). In our model, we consider two 
distinct accretion modes separated at a rate of 1% of M^dd- In the 
first state, for m ^ 0.01, we assume that the gas forms an accre- 
tion disc whose physics is adequately described by the radiatively 
efficient thin-disc model of Shakura & Sunvaevl d 1973h . The bolo- 
metric luminosity of a thin disc, Lboi, is linked to M through the 
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standard expression 



= eMc 



(5) 



When the accretion becomes substantially super-Eddington 
(Lho\ 1? yLEd d), the bolometric luminos ity is limited to 77 [1 + 
ln(rh/?7)]LEdd (Shakura & Sunvaevll 19731) , where r\ is an ad hoc 
parameter equal to 2 that allows a better modelling of bright end of 
the LF (see Section [5}. However, we do not restrict the accretion 
rate if the flow becomes super-Eddington. 

The second accretion state, with m < 0.01, is associated 
with sub-Eddington accretion flows with very low density. For such 
low accretion rates, the gas flow is unable to cool efficiently since 
radiative cooling does not balance the energy generated by vis- 
cosity. Thus, the viscous energy is trapped in the gas as entropy 
and ultimately adverted into the hole. This type of accreti on is 
known as an advection dominated accretion flow (ADAF, iReej 
ll982l : lNaravan & Yi|[l994l ; lAbramowicz et alj|l995l) . ADAFs have 
a number of distinct properties, some of which will be essential 
for the analysis in later sections (see Sections 5 and 6). For ex- 
ample, for an ADAF around a BH, only a fraction of the standard 
accretion luminosity, L = eAfc 2 , is emitted as radiation. The re- 
mainder of the viscously dissipated energy is stored in the gas as 
entropy, resulting in h ot flows w i th alm ost virial temperatures. We 
note that, as shown bv llchimarul il977l) . the ions and the electrons 
in an ADAF are not thermally coupled and, thus, reach different 
temperatures. This two-temperature virialized plasma flow is opti- 
cally thin and, for high viscosity parameters (a ~ 0.1 — 0.3), it 
acquires a quasi-spherical geometry around the BH, which resem- 
bles spherical Bondi accretion. However, the accretion is entirely 
due to dissipation via viscous forces rather than gravity. The bolo- 
metric luminosity of the flow in this case is equal to the luminosity 
emitted by the various cooling processes. For m < 10 _3 a< 2 it is 
only due to Comptonisation, whereas for m > 10 _3 o: 2 the cooling 
is split between Compton and synchrotron emission. 

As the accretion rate increases, the emission of the energy pro- 
duced by the viscous processes in the gas becomes more efficient. 
Above some critical accretion rate, m c , the radiative efficiency of 
the gas is so high that the flow cools down to a thin disc. The crit- 
ical accretion is independent of the BH mass but depends strongly 
on the viscosity, m c ~ 1.3a 2 . Taking m c = 0.01 implicitly fixes 
the value of a to be 0.087 for all the ADAFs in our model. The 
lu minosity of the flow for the two different regimes is then given 
bv lMahadevanl h997h . 

/ eMc 2 [OAm/3/a 2 ] , rh > 7.5 x 10" " 



Jbol.ADAF — 



{eMc 2 [4x 10~ 4 (l-/3)] , m< 7.5x10" 



(6) 



where /3 is related to the Shakura-Su nyaev viscosity parameter a 
throu gh the relation a » 0.55(1 — 0) JHawlev. Gammie. & Balbusl 
Il995t) . The expressions in the square brackets in Eq. [6] shows how 
much less efficient the cooling is in an ADAF compared to the stan- 
dard efficiency e of a thin disc. For example, at m = 0.01 an ADAF 
is characterized by an accretion efficiency of 0.44e, less than ap- 
proximately half as luminous as a thin disc. 

4.2 BH spins and accretion efficiencies 

In lFanidakis et alj J2010l) we describe in detail how we model the 
evolution of BH spin. In brief, the spins of BHs change during 
gas accretion (whenever a starburst is triggered by disc instabili- 
ties or galaxy mergers in the starburst mode or during the cooling 
of gas from the hydrostatic halo in the hot-halo mode) and merg- 
ers with other BHs. Fig.|6]shows the evolution with redshift of the 
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Figure 6. Top: The median BH spin weighted by BH mass as a function 
of redshift for BHs in our model (solid black line). Bottom: The median 
BH spin weighted by disc luminosity for all AGN accreting in the thin-disc 
regime (solid black line). The values of the different percentiles (shaded 
regions) are indicated in both cases by the colour bars. 



median of the BH spin distribution at a given redshift for two differ- 
ent cases. Firstly, we show the median spin weighted by BH mass 
(upper panel) for all the BHs with Mbh > 10 6 Mq. The median 
shows an approximately constant trend from z — 6 to z — 2 and 
has a well defined value of ~ 0.25. In this redshift range BHs grow 
predominantly during disc instabilities (see Fig. [2j and therefore 
the evolution o f BH spin is governed by accretion. As shown by 
Fanid akis et all accretion of gas results in low spins with a typ- 
ical value of a ~ 0.2 (under the assumption that the gas is fed 
chaotically onto the BH). Hence, as indicated also by the different 
percentiles in the plot (shaded regions), the bulk of BHs acquire 
low spins. 

Below z = 2 BH mergers start to become an important chan- 
nel for growing the BH mass, especially when they are between 
BHs of similar mass, as expected following a major galaxy merger. 
However, this is a characteristic only of the most massive BHs in 
our model (Mbh > 10 9 M ). Mer gers between those BHs tend to 
increase the spin of the final remnant to v a lues a > 0.7 (see als o 
iBaker et alj|2007l: iBerti & Volonteril 120081 ; iFanidakis et al.ll2O10h . 
Eventually, at low redshifts BH mergers give rise to a population 
of rapidly rotating BHs. The appearance of these BHs, as indicated 
also by the different percentiles below z — 2, increases the dynam- 
ical range of the predicted spins up to values of 0.998 (the range of 
spin values in the distribution of BHs in Fig. [6] is smaller because 
we show only up to the 90 th percentile of the data). The environ- 
mental dependence of BH spins is illustrated in Fig.|7] where we see 
that the most rapidly-rotating BHs populate the centres of the most 
massive DM haloes. In contrast, slowly rotating BHs are found in 
low mass halo environments. This is a consequence of the correla- 
tion between BH mass and spin in our model: rapidly rotating BHs 
have masses > 10 9 Mq. These BHs are hosted by massive ellip- 
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Figure 7. The distribution of rotating BHs at z = in the Millennium simulation. The two panels show the same DM distribution in a cube of comoving 
length 100 Mpc h _ 1 colour coded according to density (black represents the peaks of DM density). Over plotted with blue spheres in the left plot are galaxies 
with rapidly rotating BHs (a > 0.7) predicted by the semi-analytic model. Similarly, in the right plot are shown a sample of galaxies with slowly rotating 
BHs (a < 0.7) randomly chosen from our data and equal in number to the a > 0.7 BHs. The size of the spheres is proportional to the spin of the BH that the 
galaxy hosts. 



tical galaxies in our model which inhabit the centres of the most 
massive haloes. 

In contrast, the median spin of actively growing BHs in our 
model does not reveal the presence of rapidly spinning BHs. This 
is demonstrated in the lower panel of Fig.[6]where we have selected 
the actively growing BHs in our sample (m > 0.01) and plot the 
median of their spin distribution weighted by disc luminosity. As il- 
lustrated in the plot, the predictions are consistent with low spins at 
all redshifts even at z ~ 0. This is because in this sample we have 
selected only the actively growing BHs; we find that almost exclu- 
sively these are < 10 Mq BH and thus have low spins. Hence, 
these BHs dominate the sample and determine the over all trend of 
the median. 

The a — Mbh correlation can be further elucidated through 
the dependence of the accretion efficiency on the mass of actively 
growing BHs. This is illustrated in Fig. [8] where we show the me- 
dian of the e distribution for the sample of actively growing BHs 
(solid lines), again weighted by the disc luminosity. At z = 0, the 
median is approximately constant for Mbh 10 8 , with a typical 
value of ~ 0.07. As implied by the 10 — 90 percentiles (dotted 
lines), the dynamical range of the typical e values is very small and 
restricted to the range 0.06 - 0.08 (a = 0.1 - 0.5). For higher BH 
masses the efficiency can reach significantly higher values. This is 
a manifestation of the fact that high mass BHs have high spins and 
thus high accretion efficiencies when they accrete in the thin-disc 
regime. Fig. [8] also demonstrates that the dependence of the effi- 
ciency on the BH mass does not change with redshift. Hence, at 
all redshifts BHs have very well determined accretion efficiencies. 



It is therefore, only the accretion rate that regulates the luminosity 
output from an accreting BH. 

Note that, BHs with higher efficiencies (e > 0.15) are still 
found in our BH populations. However, they are usually very mas- 
sive (A/bh > 5 x 10 8 Mq) and do not undergo significant quasar 
activity. These BHs can support the formation of very strong jets in 
the presence of an ADAF and establish the host galaxy as a radio- 
loud A GN (high spins and B H mass are essential for high jet lumi- 
nosities; Fanidakis et al. 2010). In this case, a plot similar to the one 
in the top panel of Fig[6] weighted by jet instead of disc luminosity, 
will unveil a significant population of AGN with rapidly rotating 
BHs. Hence, BH spins (and efficiencies) can display different dis- 
tributions depending on the AGN population we are probing. 

4.3 The distribution of the AEdd parameter 

Having explored the evolution of the BH mass and spin, and cal- 
culated the accretion efficiencies for the BHs accreting in the thin- 
disc regime, we now investigate the disc luminosities of the accret- 
ing BHs predicted by the model. We calculate the bolometric disc 
luminosity in the ADAF and thin-disc regime using the formula- 
tion described in Section |4~T1 It is useful to scale Lbot, in units of 
£<Edd in order to remove the dependence of the luminosity on the 
BH mass. For this purpose we introduce the Eddington parameter, 
AEdd, defined as, 

{7m 2 , ifm<0.01 
rii, if m ^0.01 

ri[l + ]n(m/rj)], if Lboi > ^Edd 

(7) 
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Figure 8. The median of the accretion efficiency, e, as a function of BH 
mass at different redshifts (solid lines) for BHs that accrete in the thin disc 
regime. Also shown are the 10 — 90 percentiles of the e — Mbh distri- 
bution. We also indicate the values of e that correspond to spin values of 
0, 0.2, 0.4, 0.6 and 0.8 (dotted grey lines). The model predicts that low 
mass BHs have well defined low accretion efficiencies (e ~ 0.06 — 0.08) 
at all redshifts. In contrast, in the low redshift universe, a population of very 
massive BHs with high radiative accretion efficiencies (e > 0.08) emerges. 



where 7 = 0.4/3 /a 2 . 

In Fig.[9]we plot the distribution function of AEdd at different 
redshifts for all accreting sources with Mbh > 10 6 Mq. The plane 
of the plot is divided into two distinct regions: the thin-disc regime 
and the ADAF regime (note the discontinuity due to the different 
dependence of AEdd on m in the ADAF and thin disc regime). This 
distinction is essential since it will help us to readily unravel the 
space density and evolution of luminous and under-luminous AGN 
in our model. 

The first important property that is unambiguously depicted 
by Fig. [9] is the bimodal nature of the Asdd distribution. The bi- 
modality is exhibited nearly at all redshifts. The first peak falls in 
the ADAF regime and it shifts from log AEdd — —2.5 at 2 ~ 6 
to log Asdd — —4 at z ~ 0. The objects that contribute to that 
mode are all the AGN accreting during the hot-halo mode and those 
AGN accreting during the starburst mode with logm < —2. The 
second peak is located in the thin-disc regime and in this mode we 
find the most luminous objects in our model. Their peak space den- 
sity shifts from logAEdd — —0.5 at z ~ 6 to logAEdd — —1 
at z ~ 0, where it nearly disappears. The objects populating the 
second mode are AGN accreting exclusively in the starburst mode. 
Since the thin-disc regime is radiatively efficient these objects are 
expected to dominate the LF of AGN in all bands (except the radio 
and perhaps the hard X-rays). 

The relative space density of the objects in each mode changes 
dramatically with redshift. At high redshifts, AGN in the thin-disc 
regime are much more numerous than those in the ADAF regime. 
The BHs in these AGN have low masses, accrete near the Edding- 
ton limit and double their mass several times within a few Gyrs. In 
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Figure 9. The distribution of AEdd = ^bol / ^Edd for different redshifts as 
indicated by the colour bar. The plane is divided into the radiative-inefficient 
regime of ADAFs (shaded region) and radiative-efficient regime of thin 
discs. The denser shading denotes the region of the plane where the dis- 
continues transition from the ADAF to the thin-disc regime takes place. 



contrast, BHs in the ADAF regime are more than an order of mag- 
nitude less numerous. Their number density is dominated by BHs 
that have reached by that time a high mass (10 8 - 10 9 M Q ) and 
their host galaxy undergoes a disc instability or minor merger event 
that provides gas with a low accretion rate. The first BHs accreting 
during the hot-halo mode also contribute to the objects populating 
the ADAF regime. 

This picture changes in the low redshift universe. As redshift 
decreases more haloes enter the quasi-hydrostatic cooling regime 
and thus more BHs start to accrete in the ADAF regime (see also 
the evolution of the hot-halo accretion channel in Fig.f2). In this 
mode we also find AGN powered by gas accretion during gas-poor 
disc instabilities and galaxy mergers. Eventually, at z = BHs ac- 
creting via an ADAF are much more common than those accreting 
via a thin disc. 



4.4 The L h 



Mbh correlation 



When examining the distribution function of Asdd in Fig. [9] we 
find a small population of super-Eddington AGN present at all red- 
shifts. The number density of these AGN drops sharply with in- 
creasing Asdd due to the logarithmic dependence on the accretion 
rate. Whether these objects represent the most luminous AGN at 
some given redshift is not obvious since Asdd is independent of 
the BH mass. To unravel the relation between the BH mass and ac- 
creted luminosity we plot in Fig.[l0]the median of the Lboi — Mbh 
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distribution and its associated percentiles at z = 0.5, 1 and 2. 
We show predictions both for objects accreting in the thin-disc 
(lower panels) and ADAF (upper panels) regime. To guide the 
reader we also distinguish with different shading the region where 

Lboi ^ iEdd(MBH). 

When we consider BHs accreting in the ADAF regime we find 
that the Lboi — Mbh correlation is characterized by a floor at lu- 
minosities < 10 42 erg s~\ In this regime, BHs have AEdd — 1CP 4 
and therefore they comprise the majority of accreting BHs in our 
sample (remember the location of the peak of the AEdd distribution 
function in Fig.[9}. For example, at z = we find that ~ 93% of the 
BHs produce luminosities fainter than 10 42 erg s _1 . For brighter 
luminosities, we find a strong correlation between Mbh and Lboi 
which indicates that the most massive BHs must have higher ac- 
cretion rates compared to the lower mass ones. Indeed, when iden- 
tifying the accretion properties of the most massive BHs we find 
that the > 10 8 M BHs have A Ed d > 10~ 3 . Thus, these BHs are 
able to produce a significant luminosity even though they accrete 
in the ADAF regime. In fact, in our sample we find ~ 10 Mq 
BHs accreting at rn ~ 0.01 which produce luminosities as high 
as 10 46 ergs _1 . However, these BHs are very rare: at z — 0.5, 
where the space density of these BHs peaks, we find only a hand- 
ful of them (space densities lower than 10 -8 Mpc -3 ). For higher 
redshifts, their space density declines and as a consequence the 
maximum disc luminosity produced in the ADAF regime is also 
reduced. 

In the thin-disc regime, we find that accreting BHs typically 
produce luminosities greater than 10 42 ergs~ 1 . The log Lboi — 
log Mbh correlation in this regime increases monotonically until 
the slope becomes significantly shallower near the highest lumi- 
nosities achieved at a given redshift. The nature of the break in 
the slope is determined by the AGN feedback prescription in our 
model. When massive haloes reach quasi-hydrostatic equilibrium 
they become subject to AGN feedback that suppresses the cooling 
flows; some of the mass which would have been involved in the 
cooling flow is instead accreted onto the BH. In these haloes we 
find the most massive BHs in our model (> 10 9 Mq, see Fig[5p). 
Therefore, these BHs are expected to live in gas poor environments 
and when they accrete gas they usually do so via an ADAF disc. In 
this regime, the suppression of cooling flows forces the Lboi — Mbh 
correlation to evolve only along the Lboi axis since accretion via a 
thin disc onto > 10 9 Mq BHs becomes very rare. 

The correlation between Lboi and Mbh found at z = 0.5 re- 
mains approximately the same at higher redshifts. Only the ranges 
of BH mass corresponding to the bulk of the AGN shift modestly 
to lower mass. This is due to the fact that in a hierarchical universe 
accretion shifts to lower BH masses at higher redshifts (Section[3). 
In addition, the break in the slope at high luminosities becomes less 
prominent since fewer haloes are in quasi-hydrostatic equilibrium 
at high redshifts. 

The most luminous AGN (> 10 46 erg s _1 ) are exclusively 
powered by super-Eddington accretion onto ~ 10 8 — 10 9 Mq 
BHs. This implies that the most luminous quasars are expected to 
be found in ~ 10 12 - 10 13 M Q DM haloes and not in the most 
massive ones (remember the Mijalo — Mbh in Fig.|5p). This is in 
good agreement with the typical DM halo mass quasars are inferred 
to inhabit as suggested by several observational clus tering analyses 
dda Angela et alll2008fclshen et alJl2009l ; lRoss et alj|2009n . In con- 
trast, accretion characterized by lower luminosities spans the whole 
range of BH masses (10 6 - 10 10 M ). 



5 THE EVOLUTION OF THE AGN LUMINOSITY 
FUNCTIONS 

5.1 Bolometric corrections and obscuration 

The LF is calculated in redshift ranges which are determined by 
the observations we are comparing with. The contribution of each 
AGN to the LF in a range z = Z\ — 22 is weighted by a factor, 



WBH = tactive/At zl! ;z 



(8) 



where At zl , Z2 is the time interval delineated by the redshift range 
21 — Z2 and factivc is the time during which the AGN is "on" in 
the At zl]Z2 interval (in principle tactive = iacc only if the whole 
period of accretion falls within z\ — z-z). The bands for which 
we present predictions are the B-band (4400A), soft X-ray (SX: 
0.5 - 2 keV) and hard X-ray (HX: 2-10 keV). The bolometric 
corrections considered here are approximated by the following 3 rd 
degree polynomial relations (Marc oni et al] 20041) . 



log(L H x/L bo i) = 
log(Lsx/L bo i) = ■ 
log(^ B L„ B /Lboi) = 



-1.65 - 0.22£ - 0.012£ 2 + 0.0015£ 3 



1.54 - 0.24£ - 0.012£ 2 



0.0015£" 



-0.80 + 0.067£ - 0.017ZT + 0.0023£ d 



(9) 

where £ = log(Lboi/L©) — 12. lMarconi et al J derive these correc- 
tion based on a spectral template where the spectrum at A > lfim is 
truncated in order to remove the IR bump. In this way, the spectrum 
corresponds to the intrinsic luminosity of the AGN (optical, UV and 
X-ray emission from the disc and hot corona). We apply Eqs.[9]to 
both thin discs and ADAFs, though we note that there is evidence 
for a change in these corrections with AEdd dVasudevan & F abian 
120071) . 

There is evidence that the fract ion of obscured AGN decreases 
with i ncreasing X-ray luminosity dUeda et al. 2003; Steffe net al.l 
120031; IHasingeil |2004 lLa Franca et all l2005h a trend found also 
by ISimpsonl ( 120051) in a sample of broad and narrow-line AGN 
from the SDSS. The question of whether the fraction of obscured 
AGN depends also on redshift is more uncertain. If gas or dust 
in galaxies provide the obscuring medium then its abundance 
should evolve in a fashion similar to the SFR history. This will 
result in a fraction of obscured AGN that is redshift dependent. 
In addition, a strongly evolving population of obscured AGN 
is required by AGN population synthesis models to reproduce 



Gilli, Risaliti, & Salvati 1999; Ballantvne, Everett, & Murray 


2006|; Ballantvne et all 


2006|; iGilli. Comastri. & Hasingerl |2007|). 


UedaetalJ d2003) and 


Steffenet al. (2003) suggest that such a 



trend is not cle ar in AGN sampl e s of d eep X-ray surveys. However, 
the analysis bv lTreister & Urrvl d2006h on AGN samples of higher 
optical spectroscopic completeness indicates that the relative 
fraction of obscured AGN doe s incre ase with redshift. 

More recently, Hasinger] d2008h showed, based on a sample 
of X-ray selected AGN from ten independent samples with high 
redshift completeness, that the fraction of obscured AGN increases 
strongly w ith decrea sing luminosity and increasing redshift. Ac- 
cording to IHasingeil the dependence of the obscured fraction of 
AGN, /obsc, on Lhx, can be approximated by a relation of the 
form, 



fohS' 



-0.281 log(L H x) + A(z) 



(10) 

46 erg s -1 .A 



where A(z) = 0.308(l+2)°' 48 andL HX = 10 42 -10 
broken power law can also be used to describe A(z) which provides 
a fit with modestly better statistical significance. The best fit gives 
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Figure 10. The median of the Mbh — ^bol distribution at z = 0.5, 1 and 2 (black solid lines) for BHs accreting via a thin disc (lower panels) and ADAF 
(upper panels). The different percentiles of the distribution are colour-coded according to the bar on the right. The shaded regions represent the regime where 
the accretion becomes super-Eddington. 



a power law of the form A(z) oc (1 + z) that saturates at z = 
2.06 and remains constant thereafter. 

Hence, in order to determine the correct population of AGN in 
a luminosity bin of a given band we need to take into account the 
effects of obscuration. To do so, we utilise the d ependence of the 
obscured fraction on the Lhx luminosity found in lHasingej [2008) 
as follows. We calculate the fraction of visible AGN, / v i s = 1 — 
/obsc, at a given luminosity bin in the 2 — 10 keV band using Eq.UOl 
and then we associate the value of / v i s with the corresponding B- 
band or soft X-ray luminosity bin using Eqs.[9] The LF can then be 
expressed as, 



dlog(Lx) 



/vis(iHX, z) 



dlog(L;> 



(11) 



In our analysis, we choose the single instead of the broken power 
law fo rm of A(z) to allow for a comparison with previous work 
(e.g.. lLa Franca et alj|2005h . In this way we provide a simple, yet 
well constrained, prescription for the effects of obscuration in the 
AGN of our model. 



5.2 The optical LF 

In Fig.Qj]we present the bj-band LF of quasars in 9 different red- 
shift bins between 0.4 < z < 4.25. Our predictions are shown 



before and after applying the obscuration effect (dashed and solid 
black lines respectively). Absolute magnitudes are first calculated 
in the B band using the standard expression 



Mb = -10.44 - 2.51og(L B /10 40 erg s" 



(12) 



for magnitudes in the Vega system and then conver ted to the bj 
band u sing the correction M^ } = Mb — 0.072 dCroom et al.l 
2009b). Predictions are shown for the FanlOb and for comparison , 
the predictions from the modified version of the Bower et al. ( 200i 
model (solid orange lines) p resented in iFanidakis et alj J201 
Fan 10a). The iFanidakis et alj model was constrained to match the 
observed LF of quasars at z ~ 0.5. In brief, the model assumes a 
constant obscuration fraction of /„i s = 0.6 and a bj-band bolomet- 
ric correction of 0.2. Finally, our predictions for the 2 = LF, af- 
ter applying the obscuration effect, are shown in every redshift bin 
(black dotted lines) to help us assess the evolution with redshift. 

Our predictions are compared to the 2SLAQ+SDSS QSO LFs 
JCroom et alT20 09b). QSO magnitudes in the 2SLAQ sample were 
obtained in the g-band and therefore need to be conve rted to the bj 
band c onsidered here. We use Mbj = M g + 0.455 dCroom et al.1 
2009b) where the K-corrections for the g band have been normal- 
ized to z — 2. In addition, we plot the 2QZ LFs bv lCroom et al.l 
d2004l) . The observed LFs agree very well with each other particu- 
larly at bright magnitudes (Mb , < —24). The modest disagreement 
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Figure 11. The 6j-band quasar LF in the redshift range 0.4 < z < 4.25. Predictions are shown for the model described in this paper (FanlOb; solid black 
lines) and the FanlOa (solid orange lines) model. For reference, we show the 2 = LF in the FanlOb model (dotted black line). We also show the evolution 
of the LF in the FanlO b model befor e applying the obscuration (dashe d black line s). The observations represent the LFs estimated from the 2SLAQ+SDSS 
(blue-filled circles, Croom et al. 2009b), 2QZ (red-filled circles, ICroom et alj2004l) and SDSS (brown-filled triangles, Richards et al. 2006) surveys. 
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Table 3. Typical values of / v ; s in three magnitude bins and its evolution 
with redshift. 



M bj = -20 M bj = -24 M bj = -28 



z 


= 0.4 


53.3% 


84.3% 


100% 


z 


= 0.68 


49% 


81% 


100% 


z 


= 1.06 


44.9% 


76.9% 


100% 


z 


= 1.44 


42.2% 


73.3% 


100% 


z 


= 1.82 


27.8% 


69.9% 


100% 


z 


= 2.20 


24.7% 


66.7% 


98.7% 



seen at the faintest magnitudes is due to the different K -correction 
applied to the 2QZ QSO sample bv lCroom et al.l ( 12004) . For higher 
redshifts, we inc lude the LF derived f rom the SDSS Data Release 
3 sample (DR3. iRichards et alJ^OOfj) . The SDSS DR3 LF is ob- 
tained in the i band, therefore we need to convert it to the bj band. 
WeuseMj,, = Mj + 0.71, assuming a spectral index of a v — —0.5 
dCroom et ai]|2009bh . 

The AGN model contains one free adjustable parameter that 
must be constrained by observational data. This is the / q parameter 
in Eq. [2] which sets the accretion timescale. We choose the fej band 
for constraining the value of / q because it is the most sensitive to 
the modelling of AGN. Other free parameters are constrained by 
the galaxy formatio n model and are adjusted to reproduc e the ob- 
served galaxy LFs teower et al.ll200tj : lLagos et alJl2O10L see also 
iBower et afl [2010! for further discussion of our parameter fitting 
philosophy). Since / q is essentially the only free parameter in our 
model, the predictions presented in these Sections are genuine pre- 
dictions of the underlying galaxy formation model. By fixing / q to 
10 we obtain an excellent overall match to the observed QSO LFs 
in the 0.4 < z < 2.6 redshift interval. For higher redshifts, our 
model provides a good match, however, it predicts a much steeper 
slope for the bright end compared to the SDSS DR3 LF. 

A comparison between the predictions for the z = 0.4 — 4.25 
and z — LFs in Fig. QT] shows that quasars undergo significant 
cosmic evolution. For example, quasars with M bj = — 25 increase 
in space density from ~ 10~ 8 to ~ 10~ 6 Mpc -3 between z — 
and 2 = 2.2 — 2.6. This strong evolution in the space density of 
quasars is due to the fact that disc instabilities and galaxy merg- 
ers, the two processes that trigger accretion during the starburst 
mode, become more frequent at higher redshifts (see Fig. [TJ. Yet, 
the strong evolution of quasars does not affect only their space den- 
sity. At a fixed space density of 10~ 8 Mpc -3 the quasar LF bright- 
ens from M bj ~ -25 at z = to M b , ~ -28 at z = 2.2 - 2.6. 
This is attributed to the fact that the cold gas becomes more abun- 
dant with increasing redshift. This is equivalent to a strong increase 
in the gas reservoir available for feeding the central BHs since more 
gas is turned into stars when a starburst is triggered. 

The processes that are responsible for driving the formation 
of stars in starbursts are also responsible for the cosmic evolution 
of quasars. Qualitatively, the strong link between the formation of 
stars and quasar activity (and therefore BH growth) can be illus- 
trated by considering the FanlOa model. In Fig[T]we can see that the 
SFR density in bursts in the FanlOa model shows a milder evolu- 
tion with redshift compared to the FanlOb model. Hence, the model 
predicts considerably that less gas is available for accretion which 
then results in more modest evolution of the quasar LF. As a con- 
sequence, the predictions of the model FanlOa provide an overall 
poor match to the observed LFs. 

The evolution of the SFR density with redshift does not im- 



ply a proportional evolution of quasar luminosities on the M bj — 
$(Mbj) plane. The space density of the brightest quasars increase 
with redshift. However, the density of objects around the break 
in the LF increases more quickly, leading to a steepening of the 
LF as illustrated by FigQT] This differential evolution is deter- 
mined exclusively by the accretion physics in our model. In the 
0.01 < rh < 1 regime, the disc luminosity scales in propor- 
tion to rh. However, when the flow becomes substantially super- 
Eddington, the luminosity instead grows as ln(rh) and therefore 
the dynamical range of predicted luminosities decreases dramati- 
cally. The logarithmic dependence of luminosity on the accretion 
rate has a strong impact on the shape of the LF, resulting in a very 
steep slope at the bright end. This becomes apparent in the highest 
redshift intervals (z > 1.44) where the relative number of super- 
Eddington accreting sources becomes significantly higher than in 
the lower redshift intervals (because more gas is available for ac- 
cretion, see Fig. |9). Since these sources exclusively populate the 
brightest magnitude bins, the bright end of the LF now becomes 
significantly steeper. 

Another factor that influences the evolution of the quasar LF is 
obscuration. Low luminosity quasars are heavily obscured accord- 
ing to our obscuration prescription, and thus remain well buried 
in their host galactic nuclei. This is clear from Table [3] where we 
summarize the evolution of the / v i s value in the redshift range of 
interest. In the highest redshift intervals the obscuration becomes 
more prominent affecting even the brightest sources. However, the 
intense accretion activity during the starburst mode dominates the 
evolution of the faint end, and therefore a strong cosmic evolution 
in the space density of the faintest sources, similar to that of the 
brightest sources, is still observed. Nonetheless, the obscuration 
significantly influences the cosmic abundance of the quasar popu- 
lations dominating the faint end, something that will become more 
evident in Section [531 where we will explore the cosmic evolution 
of quasars of given intrinsic luminosity. 

Finally, we point out that those AGN that are powered by 
an ADAF contribute significantly to the space density of AGN 
fainter than M bj ~ —22 in the low redshift universe. In fact, 
in the range 0.4 < z < 0.68, we find ADAF sources with 
rh ~ 0.001 — 0.01 contributing also to the knee of the LF. These 
systems represent the rare > 10 9 Mq accreting BHs in Fig. [10] 
Identifying the ADAF systems with optically bright quasars intro- 
duces a caveat for the model since quasars typically show high ex- 
citation spect ra which indicate the presence of a b right, UV thin 
disc (see e.g.. lMarchesini. Celotti. & Ferraresd2004T) . Nonetheless, 
these systems accrete where the transition to a thin disc takes place. 
Observations of stellar mass BH binary systems show that this tran- 
sition is complex, probably taking on a composite structure with 
the thin disc replacing the hot flow at progressively sma ller radii 
(see e.g,. the review by iDone. Gierlinski. & Kubota 2007). Such a 
configuration could possibly produce high excitation lines while re- 
taining also the ADAF characteristics. For a more comprehensible 
representation of the objects populating the LF, we refer the reader 
to Fig.[l4]in Section [54l where we decompose the LF into the con- 
tribution from ADAF, thin-disc, and also super-Eddington sources. 

5.3 The X-ray LFs 

The good agreement between our model predictions for the optical 
LF and the observations motivates us to study also the evolution 
of X-ray AGN as a function of redshift and intrinsic luminosity. 
X-rays account for a considerable fraction of the bolometric lu- 
minosity of the accretion disc and therefore provide an ideal band 
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Figure 12. The predictions of our model for evolution of the soft X-ray LF (solid black lines). Predictions are shown before (total: dashed lines) and after 
(visible: solid lines) we apply the e ffects of obscuratio n. In addition, we quote the total fraction of AGN with Lgx > 10 42 erg s _1 that are visible in every 
luminosity bin. Data are taken from Ueda et alj fe003l) . 
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Figure 13. The evolution of the hard X-r ay LF in our model (solid black line s). Predictions are shown only for z < 1.6 and no obscuration correction is 
applied. Observational data are taken from Hasin ger, Mivaii. & Schmidt (2 0051) . 
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for studying the cosmic evolution of AGN. The observed contin- 
uum spectra of AGN in the X-rays can be represented by a simple 
power law in the hard part of the spectrum (2 — 10 keV) with 
a universal spectral index of ax ~ 0.7 which is independent of 
the luminosity over several decades. In the soft part of the spec- 
trum (0.5 — 2 keV) the power law continuum is strongly attenuated 
by absorption which is typically in excess of the Galactic value 
and corresponds to a mediu m with hydrogen column densities of 
J Vh g 10 20 - 10 23 cm" 1 teevnoldd[l997l ; iGeorge et alj|l998l ; 



IPiconcelli et~aT]|2005r) . 

The cosmic evolution of X-ray AGN has been in- 
vestigated by employing AGN selected both in the soft 
(Mivaii, Hasinger, & Schmidt 2000; Hasinger, Mivaii .~& Schmidtl 
2005b and hard p a rt of the spectr u m foeda et alj 120031 ; 



La Franca et ail |2005t iBarger et all 120051 ; lAird et al J |201U> . Mir 



roring the strong evolution of the optical LF, the LF of X-ray se- 
lected AGN also evolves strongly with cosmic time. The LFs in 
soft and hard X-rays have similar characteristics, evolving at the 
sa me rate but differ significantly in normalization. For example, 
the|Hasinger, Mivaii, & Schmidt (2005]) soft X-ray sample has ap- 
proximately a 5 times smaller global LF normalization com pared 
to that of the hard X-ray LF estimated bv lUeda et al] ( l2003h . This 
difference is most likely attributed to the fact that a large fraction of 
AGN is obscured in the soft X-rays (ty pe-2 AGN are by a factor of 
4 mo re numerous than type-1 AGN, Risaliti, M aiolino. & Salvatil 
1 19991) . 

We present in this section our predictions for the X-ray 
LFs. We calculate the X-ray emission in the 0.5 — 2 keV and 
2—10 keV bands using Eq. {5}. Our predictions for the soft 
and hard X-ray LFs are shown in Figs. [12] and [T3] respectively 
(so lid lines). The predictions are comp ared to the LFs estimated 



by Hasinger, Mivaii, & Schmidt (2005) for the soft X-rays and 
lUeda etaljfeool) for the hard X-rays. Note that we do not show es- 
timates fro m observations that consist on ly of upper limits (empty 
bins). The lHasinger. Mivaii. & Schmidtl sample comprises unab- 
sorbed AGN and therefore we need to take this into account for 
the samp le in the soft X-r ays the effect of obscuration. We do so by 
using the Hasingei] J2008h prescription as explained in the previous 
section. To quantify the effect of absorption we also plot in the soft 
X-rays the total population of AGN (dashed lines). 

Overall, our model provides a very good match to the observa- 
tions in both soft and hard X-rays. We find discrepancies between 
the observations and the model predictions in the redshift range 
z — 0.4 — 0.8, where the space density of the < 10 43 erg s _1 soft 
X-ray AGN is under predicted. A similar disparity is also seen in 
the z = 0.8 - 1.6 bin for the 10 43 - 10 45 erg s" 1 AGN. However, 
it is important to bear in mind that our predictions in the soft X- 
rays depend strongly on our prescription for the obscuration. Since 
our predictions for the total AGN population are always well above 
the observations, it is likely that these discrepancies are attributed 
to the modelling of the obscuration. We also find that the model 
predicts at high redshifts a somewhat steeper bright end compared 
to the observations, although the very brightest point could be due 
to beaming effects i n a small fraction of th e more numerous lower 
luminosity sources iGhisellini etai]|2010l) . Nonetheless, the very 
steep slope in our predictions is again a manifestation of the Ed- 
dington limit applied to the super-Eddington sources, an effect that 
becomes more significant at higher redshifts as already explained 
in the previous section. 

As illustrated by the predictions for the soft X-ray LF in 
Fig. [12] the space density of the faintest soft X-ray AGN in our 
model is strongly affected by the obscuration. As illustrated by the 



LFs in the soft X-rays, AGN with L S x ^ 10 42 - 10 44 erg s _1 are 
heavily obscured. In fact, in the redshift interval z — 0.015 — 0.2 
only 19% of the total AGN population with L S x ^ 10 42 erg s~ x 
is visible in soft X-rays. Henc e, the obscured AGN outnumber th e 
visible AGN by a factor of 4 dRisaliti. Maiolino" & Salvatii ri999). 
This fraction of visible AGN becomes even smaller at higher red- 
shift because the obscuration becomes more prominent with red- 
shift. At z > 3.2 only a negligible fraction of the AGN are seen in 
soft X-rays. Based on these results we estimate that approximately 
10% of the total number of AGN with Lsx 10 42 erg s _1 are vis- 
ible in the soft X-rays in the z = 0.015 — 4.8 universe. We note 
that these estimations are based on the obscu r ation model which 
according to the AGN sample from Hasinger (2008) is evaluated 
only for luminosities 10 42 erg s -1 < Lhx < 10 46 ergs _1 . For 
luminosities lower than 10 42 erg s _1 we assume that / v i s remains 
constant and equal to / v i s (10 42 erg s _1 ). It is therefore likely that 
we underestimate the value of / v i s - 



5.4 The bolometric LF 

We have so far studied the LF of AGN and quasars in different 
wavebands using the bolometric corrections in Eq. {9]l and com- 
pared them to the available observations. We find that our predic- 
tions match reasonably well the observations, particularly in the 
fej band. Further insight into the evolution of AGN can be gained 
by calculating the bolometric LF and studying its evolution with 
redshift. This will provide an overall characterisation of the global 
AGN population and might reveal additional trends which are per- 
haps unseen when exploring the LF in the other bands we studied 
earlier. 

To calculate the bolometric LF we consider all accreting ob- 
jects (in both the starburst and hot-halo mode) taking into account 
the regime these objects accrete in (ADAF or thin disc). This time 
we probe a wider baseline in luminosity and thus we expect to 
see clearly the contribution of the AGN powered by ADAFs. We 
calculate the bolometric LF at some redshift z following the tech- 
nique described in the previous sections. AGN are sampled over 
a period equal to 20% of the age of the universe at redshift z. 
We need to stress that the resolution of the Millennium simula- 
tion has an impact on the gas properties of the < 10 44 erg s _1 
objects at high redshifts. Therefore, we compare the N-body re- 
sults with high-resolution simulations using Monte-Carlo (MC) 
halo merger trees to test the limit of our predictions. The MC al- 
gorithm we use to generate the DM halo m erger trees has been 
presented in Parkins on. Cole. & Hellvl J2008h . The algorithm is a 
m odification of the Extended Press-Schechter algorithm described 
in Col e et al. and accurately reproduces the conditional MFs 

predicted by the Millennium N-body simulation. The comparison 
shows that the calculation with N-body trees becomes incomplete 
for Lboi < 10 44 - 10 45 erg s _1 AGN at z > 4. Therefore, at high 
redshifts we show predictions only down to the resolution limit of 
the Millennium simulation (the limit is indicated by the vertical 
dotted lines). 

Our predictions for the bolometric LF are shown in Fig. [14] 
(solid black lines) in the redshift range z — 0.2 — 6. In addi- 
tion, we show independently the contribution to the LF from the 
ADAF, thin-disc, and Lboi ^ ^Edd sources (dashed coloured 
lines). Our results ar e compared to the bolometric LF ext racted 
from observations by Hop kins. Richards. & Hernquist] J2007l) . us- 
ing a combination of a large set of observed LFs in the optical, 
X-ray, near- and mid-infrared wavebands. The authors are able 
to reproduce the bolometric LF and the individual bands by em- 
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Figure 14. The bolometric LF predicted by our model in 6 redshift bins between z = 0.2 and 2 = 6. Also shown is the bolometric LF es- 
timated by Hopk ins, Richards & HernquisJ J2007t) using a large sample of infrared (IR) optical, soft and hard X-ray LFs (we refer the reader to 
Hopkins, Richards, & Hernquis J I2007I for a description of the observational samples used by the authors). The different bands are depicted by different 
colours as indicated by the labels. The dashed lines show the contribution to the LF of the ADAF (blue), thin disc (red) and super-Eddington sources (green). 
The vertical dotted lines indicate the resolution limit of the Millennium simulation. 



ploying their best-fit estimates of the column density and spec- 
tral energy distribution, but also a prescription for the obscura- 
tion fraction which is assumed to be a function of the luminosity 
(though not of redshift). Their AGN LF also show the downsiz- 
ing of AGN activity at low redshifts, which is expressed through 
the stee pening of the faint-end slope b elow z = 1. The data 
used by Hopkin sT Richards. & Hernquistl cover redshift wide red- 
shift ranges and sometimes the same data sets appear in more than 
one bin (especially the soft X-ray data at z > 0.8). Therefore, com- 
parison between the data and our predictions should be treated with 
caution because the observations over such wide redshift bins could 
hide evolutionary trends. 

Nevertheless, our predictions for the bolometric LF match 
very well the observational estimates across a wide range of red- 
shifts. In the low-redshift universe (z < 1), the model reproduces 
the faint and bright ends of the observed LF remarkably well. At 
the faint end, the model predicts a significant contribution from 
AGN powered by ADAFs. These are predominately massive BHs 



accreting at low accretion rates during the hot-halo mode and they 
account for the increasing abundance with decreasing redshift of 
the faint AGN (see also Fig. |5J. By contrast, the bright end is 
always populated by AGN radiating near or greater than the Ed- 
dington limit. These are the AGN that harbour the rapidly growing 
10 7 - 10 8 M BHs in the low-redshift universe. 

At high redshifts (z > 1), the contribution from the ADAF 
sources becomes less important. This is because not only fewer 
haloes reach quasi-hydrostatic equilibrium, but also most of the 
accretion during the starburst mode takes place in the thin-disc 
regime. The LF at these redshifts is dominated by AGN whose BHs 
accrete during the starburst model. The BHs in these AGN grow 
very fast as they usually double their mass several times during a 
single accretion episode. This gives rise to a strongly evolving pop- 
ulation of super-Eddington sources that dominate the bright end of 
the LF. 

The bright end of the z 1 bolometric LF has a very steep 
slope compared to the observational estimate, indicating a defi- 
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Figure 15. Our predictions for the bolometric LF at z = 2 considering four 
different values of the accre tion timescale parameter, f a , as indicate d in the 
legend. Data are taken from Hopkin s, Richards. & Hern quisl (2007). 

ciency of very bright quasars (Lbdl > 10 48 ergs _1 ). The inabil- 
ity of the model to predict the luminosities of the most luminous 
quasars is due to the fact that the accretion luminosity is constrained 
by the Eddington limit. As discussed before, applying the Edding- 
ton limit to the accreting sources in our model results in a very 
steep decline in the space density of the super-Eddington sources, 
an effect clearly seen earlier in the optical LE The reader needs 
to keep in mind though that, as far as the observational estimates 
are concerned, many of the brightest luminosity bins do not include 
any objects and therefore, they constitute only upper limits. 

To achieve higher quasar luminosities we need to increase the 
accretion rates in the most massive BHs to values much higher than 
the Eddington rat^B This is necessary because even if we assume 
that all our BHs accrete at the Eddington limit the maximum lu- 
minosity we can produce is ~ 10 48 erg s _1 (when we consider, 
for example, a 10 10 Mq BH) which is still lower than the highest 
luminosities observed. Note though that, AGN feedback prevents 
BHs more massive than 10 9 Mq from accreting at or near supper- 
Eddington accretion rates in our model. 

To increase the accretion rate of a BH we can either increase 
the amount of gas that is fed into it or decrease the accretion 
timescale. The former solution corresponds to increasing the value 
of /bh (the parameter that determines the fraction of gas avail- 
able for accretion) which has already been tuned to provide a good 
match to the local BH density and ME The latter solution gives 
us the freedom to adjust the accretion rates without changing the 
BH mass properties. By decreasing the value of / q in Eq. [2] we 
can obtain higher m values and therefore boost the produced bolo- 
metric luminosities to Lboi > 10 48 erg s _1 . This is illustrated in 
Fig. Q3] where we show the bolometric LF at z = 2 assuming 
/ q = 20, 10, 1 and 0.5. However, changing the value of f q has 

1 For comparison, a BH accreting at m ~ 100 is only ~ 5 more luminous 
that a BH accreting at rh = 1. 



a strong effect on the faint end of the LF since it decreases the 
space density of the low accretion rate objects (the space density 
of the < 10 43 erg s _1 AGN remains unchanged since that part of 
the LF is dominated by AGN in the hot-halo mode in which the 
accretion timescale is derived directly from the cooling timescale 
of the gas). In addition, it reduces significantly the duty cycle of 
actively growing BHs by limiting the typical accretion timescale to 
~ 10 5 — 10 6 yr if / q ^ 1. Hence, having a constant value of / q 
for all the accreting BHs limits our ability to account for the entire 
Lboi baseline. Perhaps this suggests that a more sophisticated treat- 
ment of the accretion timescale is therefore needed. However, this 
is beyond the scope of the present analysis. 

5.5 The evolution of cosmic AGN abundances: cosmic 
downsizing? 

Our predictions for the optical, X-ray and bolometric LFs, sug- 
gest that AGN undergo important cosmic evolution. Their space 
density was significantly higher at earlier epochs, an evolution- 
ary trend which suggests that AGN activity in the past was much 
more intense. Interestingly, our predictions indicate that the effects 
of the different physical processes considered here (the two accre- 
tion modes and the obscuration prescription) must have an impor- 
tant, and indeed variable, effect on the evolution of faint and bright 
AGN. Therefore, to gain more insight into the cosmic evolution 
of AGN we study how the abundance of different luminosity pop- 
ulations of AGN evolves with redshift. This is shown in Fig. 1161 
where we show the cosmic evolution of optically selected AGN in 
1 1 magnitude bins (left panel) and soft X-ray selected AGN in 8 
luminosity bins (right panel). We have depicted separately the evo- 
lution of the total number of AGN (before applying obscuration; 
dashed lines) and visible population of AGN (after applying ob- 
scuration; solid lines) in each bin. In this figure we also show the 
sum over all magnitudes with Mt } ^ — 18 and luminosities with 
log Lsx *5 10 42 erg s _1 (black solid lines for the visible and black 
dashed lines for the total populations). 

A feature immediately evident in both wavebands is that the 
evolution of the space density of AGN is shallower for fainter 
sources. For example, in the optical waveband the increase in the 
space density of the Mjy = —21 AGN is characterized in the red- 
shift interval z = — 1 by a power law of the form (1 + z) f3 
with a very shallow slope of f3 ~ 0.5. The slope becomes gradu- 
ally steeper reaching a value of ~ 2.3 for the M&j = — 26 AGN. 
The steepening of the slope is evident in both the total and visible 
populations and it is driven primarily by the two distinct accretion 
channels in our model as follows. The faintest sources depicted 
in the plots of Fig. [16] comprise a combination of objects accret- 
ing during the starburst and hot-halo mode. The hot-halo mode is 
responsible for building up the majority of them in the low red- 
shift universe whereas the starburst mode accounts for them ex- 
clusively in the high redshift universe (cf. Figs. [9] and [10] but also 
Fig.ll4t. The slope of these populations changes very slowly since 
the transition from the hot-halo mode to starburst-mode domina- 
tion is not characterized by a strong change in the space density of 
AGN (the decrease of hot-halo sources is compensated by the in- 
crease of starburst sources). However, when we consider brighter 
populations, the contribution of the hot-halo mode becomes less 
important (or vanishes for the brightest populations). As a conse- 
quence, in the low redshift universe we find a much steeper slope 
for the bright AGN since their evolution is driven primarily by the 
starburst mode. 

A second noticeable feature is the strong reduction of the vis- 
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Figure 16. The cosmic evolution of different magnitude and luminosity classes of AGN selected in the optical Mb , -band (left panel) and soft X-rays (right 
panel) respectively. In addition, we show the cosmic evolution of the sum of all AGN with ^ —18 in the optical and log Lsx ^ 10 42 erg s" 1 in the 
soft X-rays (black lines). In all cases we show with dashed and solid lines the evolution of the space density of AGN before (total) and after (visible) applying 
the effects of obscuration. 



ible populations compared to the total AGN populations. The ef- 
fect is stronger in the faintest populations where we see a decrease 
in the space density of AGN sometimes greater than an order of 
magnitude. This has a significant effect on the redshift at which 
the space density of AGN peaks. For example, the space density 
oftheMfcj sC -18 and log L S x > 10 42 erg s" 1 AGN (dashed 
black lines) peaks at z ~ 3.5, which is the same redshift where 
the SF in bursts peaks (Fig[T]l. The same behaviour is seen also in 
the individual luminosity AGN populations in both the optical and 
soft X-rays. However, the visible population shows a different be- 
haviour. AGN activity, as observed in the soft X-rays for example, 
seems to peak at 2 ~ 2 — 2.5, much later than that of the total AGN 
population. More interestingly, when considering each luminosity 
population separately we notice that the brightest sources peak at 
2 ~ 3 whereas the faintest ones at z < 2. Similarly, although less 
obviously, in the optical the space density of the faintest quasars 
peaks at z ~ 2 whereas that of the most luminous increases mono- 
tonically until z ~ 3. This behaviour is ascribed to the fact that 
fainter populations are subject to stronger obscuration. Therefore 
their abundance reaches a maximum at lower redshifts and beyond 
that it declines faster. The strong obscuration also contributes to the 
further flattening of the slope making the evolution of these AGN 
quite modest. 

Trends similar to those seen in the cosmic evolution of AGN 
in our model have been reported i n the literature. For exam- 
ple, [Hopki ns, Richards, & Hernquistl J2007I) find in their analysis 



of the bolometric LF that faint AGN evolve much more slowly 
than bright AGN (Fig. 9 in their analysis). The modest evolution 
is evident also in the optica l and soft/hard X-rays. In addition, 
Hasinger, Miyaji, & Schmidt] J2005h report that in their sample of 
soft X-ray AGN, the redshift at which the space density of AGN 
peaks changes significantly with luminosity: the peak of AGN with 
10 42 erg s" 1 < L S x < 10 43 erg s" 1 is at z ~ 0.7 but gradu- 
ally shifts to higher redshifts, rea ching z ~ 2 f or 10 45 erg s _1 < 
Lsx < lO^ergs" 1 . Similarly, ICroom et alj J2009bl) found the 
same downsizing behaviour in the optical sample of quasars from 
2SLAQ+SDSS. The faintest quasars in their sample peak at z ~ 
0.6 - 0.8 (-21.5 > M g > -22.5) whereas the brightest sources 
(M g ^ —25.5) seem to monotonically increase in density up to 
2 ~ 2.5. Further evidence that the redshift at which the space den- 
sity of qua sar peaks is a strong function of luminosity form the 
analysis of ICroom et al.l is the inability of pure luminosity evolu- 
tion models (PLE; models where only the luminosity of objects 
changes) to provide an excellent fit to the data. Instead, luminosity 
and density evolution (LEDE; models where both the luminosity 
and density of objects change) are needed to provide the best fit to 
the data. 

These predictions pose the following key question. Why does 
AGN activity shift from high-luminosity objects in the high- z uni- 
verse to low-luminosity objects in the low-2 universe? Is this clear 
evidence against hierarchical galaxy formation models as has been 
interpreted by many authors? In our model, the downsizing of AGN 



© 0000 RAS, MNRAS 000, 000-000 



20 Fanidakis et al. 



is expressed through the shallow slope that characterizes the evolu- 
tion of faint AGN and the dependence on luminosity of the redshift 
at which the AGN density in a given luminosity range peaks. The 
former arises naturally from the interplay of t he starburst and hot - 
halo mode, just as in the galaxy population dBower et alj|2006h . 
Therefore, the SF and cooling processes in a hierarchical cosmol- 
ogy can account consistently for the different evolution with red- 
shift of the individual luminosity populations. The latter is driven 
largely by obscuration biases rather than the nature of their host 
galaxies. The analysis presented earlier in this section shows clearly 
that the cosmic evolution of the AGN space density is subject to 
strong obscuration which depends significantly on both luminos- 
ity and redshift. Therefore, AGN populations that are subject to 
obscurat ion could lead to mi sleading conclusions. For example, in 
both the ICroom et al.1 J2009bh OSO (these are only type-1 A GN 
ICroom et aljj2009al) ) and lHasinger. Mivaii, & Schmidt! J2005h soft 
X-ray AGN samples no correction for obscuration is adopted. 
Given that at low redshifts only u nabsorbed AGN are detected in 
soft X-rays dLa Franca et ai] |2005h . these samples can only trace 
the evolution of a small fraction (and also a specific class) of the 
total AGN population. 

Hence, to probe the intrinsic evolution of AGN one needs to 
correct for observational biases caused by obscuration (note also 
that incompleteness effects at low luminosities and high redshifts 
could also influen ce observational results about the evolution of 
AGN). lUeda" et alj ( 120031) investigated the cosmological evolution 
of hard X-ray AGN using a combination of surveys such as the 
HEAP 1, A SCA and Chandra. The sample of AGN constructed by 
Ueda et all which was compared to our predictions for the hard 
X-ray LF in Section 1531 was corrected using an absorption dis- 
tribution function that dep ends on both th e lum inosity and red- 
shift. In a similar analysis, lLa Franca et alj d2005h used a sample 
of AGN from the HELLAS2XMM survey to study the evolution 
of the hard X-ray LF taking into account selection effects due to 
X-ray absorption. In both these studies, a careful ex amination of 
the "cosmic evo lution" plots (Fig. 12 in lUeda et al.1 and Fig. 8 in 
lLa Franca et alj) reveals no evidence for downsizing in hard X- 
rays. 

When accounting for these effects, our model suggests that the 
cosmic interplay between the starburst and hot-halo mode results in 
complex evolution scheme for the AGN. In this scheme, low lumi- 
nosity AGN avoid the cosmic fate (namely becoming dramatically 
less numerous with time) of their bright counterparts, since accre- 
tion during the hot-halo mode provides gas for enabling modest 
AGN activity in the low-z universe. Note that, as already implied 
by our analysis, modest activity implies low accretion rates rather 
than low BH masses. The AGN in the hot-halo mode show a wide 
range of BH masses, and therefore, in our model the AGN down- 
sizing does not imply that the growth of low-mass BHs is delayed 
to low redshifts. Averaged-sized BHs accreting at low rates as an 
interpretation of the reported AGN downsizing is also supported 
by observations of X-ray selected AGN in the Chandra Deep Field 
South at * < 1 l lBabic et all2007t) . 



culation of BH spin. The fiducial model in this paper uses an im- 
proved SF law, as implemented by Lagos et al. (2010). 

Using this model we calculate the cosmic evolution of the fun- 
damental parameters that describe BHs in our model. These are the 
BH mass, Mbh, BH spin, a, and accretion rate onto the BH, rh 
(expressed in units of the the Eddington accretion rate through out 
this paper). We find that at high redshifts (z ~ 6) it is mainly the 
10 6 — 10 7 M© BHs accreting at rh ~ 0.3 that are actively growing. 
This picture changes at lower redshifts where we find that the accre- 
tion activity peaks for 10 7 — 10 s Mq BHs, accreting at m ~ 0.05. 
Throughout the evolution of these BHs their spin is kept low be- 
cause of the chaotic fashion with which gas is consumed. However, 
when these BHs grow to masses of > 5 x 10 8 Mq they acquire 
high spins as a consequence of mergers with other BHs. 

Knowledge of the values of Mbh, a and rh allows us to cal- 
culate the luminosity produced during the accretion of gas. To do 
so, we assume that accretion takes place in two distinct regimes: 
the thin-disc (radiatively efficient) and ADAF (radiatively ineffi- 
cient) regime. Objects accreting in the thin-disc regime are lumi- 
nous enough to account for the observed luminosities of AGN. Us- 
ing the luminosities predicted for each AGN (formed when the cen- 
tral BH experiences an accretion episode) and a prescription for 
taking into account the effects of obscuration {Hasinger 200j|) we 
calculate the LF of all accreting objects in the optical, soft and hard 
X-rays. We find an a very good agreement with the observations, 
particularly in the optical, by adjusting the value of the / q parame- 
ter, namely the proportionality factor that determines the accretion 
timescale in our model. 

The presence of two distinct accretion channels naturally 
causes a downsizing in the predicted AGN populations. These 
channels shape the evolution of the faint AGN with cosmic time, 
indicating that the faint end of the LF is dominated by massive BHs 
experiencing quiescent accretion. By contrast, the bright end is al- 
ways populated by AGN radiating near or greater than the Edding- 
ton limit. In addition, our model suggests that a significant fraction 
of AGN are obscured in optical and soft X-rays. In fact, we find that 
~ 90% of the total number of AGN in the z = 0.015 - 4.8 uni- 
verse are not visible in soft X rays. The implications of the obscura- 
tion are further revealed when we study the cosmological evolution 
of AGN of different luminosity populations. As demonstrated by 
Fig. low luminosity AGN populations are strongly attenuated 
by obscuration. As a result, their peak of AGN activity appears 
shifted to lower redshifts. 

This analysis explicitly demonstrates that hierarchical cosmo- 
logical models for galaxy formation and evolution are able to pro- 
vide a robust framework in which the evolution of AGN can be 
studied. The ability of our galaxy formation model to reproduce 
the observed LFs and account for the mechanisms that shape the 
evolution of the faint and bright AGN populations strengthens the 
powerful capabilities of semi-analytic modelling and shows that the 
level of AGN activity implied by AGN feedback is compatible with 
observations. In future studies, we aim to further study the statis- 
tical properties of AGN by exploring their spatial clustering and 
environmental dependence. 



6 CONCLUSIONS 

In this paper we have made predictions for the evolution of AGN 
using an extended version of the semi-analytic code GALFORM. 
GALFORM simulates t he formation and evo lution of galaxies in a 
ACDM cosmology. In lFanidakis et alj d2010h . we introduced a cal- 
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